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Abstract 

This paper presents the derivation of a kinetic-balance condition for explicitly 
correlated basis functions employed in senri-classical relativistic calculations. Such 
a condition is important to ensure variational stability in algorithms based on the 
first-quantized Dirac theory of 1/2-fermions. We demonstrate that the kinetic-balance 
condition can be obtained from the row reduction process commonly applied to solve 
systems of linear equations. The resulting form of kinetic balance establishes a rela¬ 
tion for the 4^ components of the spinor of an iV-fermion system to the non-relativistic 
limit, which is in accordance with recent developments in the held of exact decoupling 
in relativistic orbital-based many-electron theory. 
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1 Introduction 


Most of relativistic quantum chemistry and molecular physics is based on the 
(first-quantized) Dirac Hamiltonian p}j7]. However, unlike its non-relativistic 
counterpart, the Dirac Hamiltonian is not bounded from below and measures 
have to be taken in order to obtain correct lower bounds for the ground- and 
excited-state energies of bound states. Depending on whether the small compo¬ 
nents of the one-fermion basis spinors are included or eliminated (by some decou¬ 
pling approach [8]), methods are classified as four- or two-component methods. 
Four-component methods rely on the kinetic-balance condition for variational 
stability. This condition is well-defined for single fermions [SHE] and can there¬ 
fore straightforwardly be applied to orbital-based methods such as the Dirac- 
Hartree-Fock approach and electron-correlation methods based on it 0ZH3D]. 
For orbital-based theories with explicit correlation factors, recent work focused on 
four-component second-order Mpllcr-Plesset perturbation theory with positive- 
energy-states projection operators in combination with the one-electron kinetic- 
balance condition [31]. Li and co-workers have studied coalescence conditions for 
explicitly correlated four-component wave functions [32] but without addressing 
the issue of kinetic balance. 

A first solution to the full problem of kinetic balance for explicitly correlated 
trial wave functions was presented by Pestka and co-workers who have published 
a series of papers investigating the relativistic helium-like two-electron systems 
treated as a two-electron system in a central potential [331438] . Their solution 
is an infinite series of transformations of the individual components of the two- 
electron 16-spinor which is truncated in order to obtain an approximately kinet- 
ically balanced trial wave function. Unfortunately, little technical information 
is provided in Refs. [33tf38] and it remains unclear how such an approximate 
kinetic-balance condition can be extended to systems containing more than two 
fermions. 

In this work, we extend the pioneering work by Pestka et al. on He-like atoms [36] 
and present a scheme which allows us to derive an explicitly correlated kinetic- 
balance condition based on row reduction and a form similar to the row-reduced 
echelon form of the augmented matrix. We begin in section [2] with the pre¬ 
sentation of the theoretical background. In section [3] we apply our scheme to a 
two-electron system. Then, in section 3] we show that the correct non-relativistic 
limit is obtained. In section [5] we illustrate how the computational cost can be 
reduced for the 7V-fermion case by introducing systematic approximations to a 
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given order in the speed of light. Finally, in section [HI we demonstrate the varia¬ 
tional stability of explicitly correlated, kinetically balanced trial wave functions 
for the ground state of the He atom. 


2 Theoretical Background 


The relativistic description of electrons based on the many-1/2-fermion Dirac 
Hamiltonian provides us with a first-quantized, i.e., semi-classical formalism cap¬ 
turing essential aspects of special relativity for molecular matter m- 


2.1 The Relativistic Electron 


A single 1/2-fermion, such as an electron, may be described by the Dirac Hamil¬ 
tonian 3T) : ., IT) 

h D = col ■ p + (3mc 2 + V . (1) 


The matrices a = (a x , ct y , a z ) and (3 are defined by anti-commutation rela¬ 
tions. The most common choice that respects these relations is the standard 
representation of 4 x 4 matrices, 


cti = 


0 2 


0 2 


with i e {x, y, z } and {3 


I2 O2 
0 2 1 2 


( 2 ) 


where cr, denotes one of the three Pauli spin matrices and I 2 is the two-dimensional 
unit matrix, p = (p x ,PyiPz) T is the momentum operator, V is an operator for the 
interaction energy due to external potentials, m is the rest mass of the fermion, 
and c is the speed of light. 


It is convenient to introduce a block structure for the one-fermion eigenfunction 
-0(r), the 4-spinor, according to the 2x2 super-structure of the four-dimensional 
OLi and (3 matrices in standard representation, 


V>(r) 




10' 

_V> s (r)_ 


». 


(3) 


where T denotes the so-called large and ’s’ the corresponding small component. 
We refer the reader to the review by Esteban, Lewin, and Sere m and the 
book by Thaller 03 for a more detailed mathematical discussion of the Dirac 
Hamiltonian and its eigenfunctions. 
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The spectrum of the Dirac Hamiltonian features three distinct parts. The first 
part comprises the discrete bound states with energies between +mc 2 and —me 2 . 
The second part is the positive continuum ranging from +mc 2 to +oo. The last 
part of the spectrum is the negative continuum ranging from —me 2 to — oo. The 
negative continuum is a source of instabilities in variational calculations where 
the Rayleigh quotient, 


E[h D ,*l>(r)\ 


(Mr)\ h D IV’O)) 

(- 00 ) I i/)(r)) 


( 4 ) 


is minimized and (usually unwanted) negative-energy continuum solutions can be 
encountered if no precautions (such as projection onto positive-energy states) are 
taken into account. For basis-set expansion techniques, Schwarz and co-workers 
showed that the finite size of ordinary basis sets may pose difficulties HUE!, 
which is therefore sometimes called the ’finite-basis disease ’ m- 


An effective means of dealing with the problem of variational collapse is the 
kinetic-balance condition pillOU1211TB] which relates the large and the small com¬ 
ponent of the 4-spinor: 

ip s (r) « ^ xp l (r) . (5) 


The derivation of this relation is straightforward. The Dirac eigenvalue problem 

(h D - El 4 )^(r) = 0 (6) 

leads to a set of two linear equations for the two components of the 4-spinor 
in Eq. (J3J). After the energy spectrum has been shifted by —me 2 to match the 
non-relativistic energy scale, this system of equations reads 

(V — E)ij> l (r) + ccr • p xjj s {r ) = 0 , (7) 

(V — E — 2 mc 2 )x^ s (r) + ccr ■ p = 0 , ( 8 ) 


where the four-dimensional operator V was assumed to be a diagonal matrix with 
the same element V as diagonal entries. We only need one of the two equations 
to relate the small component to the large one. Since a ■ p has no multiplicative 
inverse, it is more convenient to choose the second equation in order to obtain an 
expression for xp s (r). After rearranging the terms, we obtain the exact relation 
for Eq. (J 8 ]) 


V>'(r) 


ccr ■ p 

(E — V + 2 me 2 ) 




( 9 ) 
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This relation depends on the energy of the system which is not known a priori 
but is one of the desired results of the problem. Eq. ([9]) can therefore not be 
applied to our problem. Now, (E — V ) is considered small compared to me 2 so 
that we may introduce the approximation 

E — V + 2 me 2 ~ 2 me 2 (10) 

to eliminate (E — V) and arrive at the kinetic-balance condition in Eq. (]5j). 
We note that this approximation step turned out to be unimportant for the 
construction of variationally stable basis-set expansion techniques applied in four- 
component orbital-based theories pan], which we assume to remain valid for the 
IV-particle theory to be developed in this work. 

Basis-set expansions which obey Eq. (J5]) provide a variationally stable parametriza- 
tion of a trial wave function for a single fermion. Eq. ([5]) may therefore be 
formulated in terms of the one-fermion model spaces [1111361,47J 

|0 G H l and \s) E (<r • p)U l C U s . (11) 


This one-fermion kinetic-balance condition can be imposed by a transformation 

0 , 


r/W _ 
‘■'KB — 


1 1 

o 

to 

i 

O' 

O 


I 2 

l- 

CM 

O 

i_ 


0 2 

(7 ■ p 

p 


( 12 ) 


(with p = |p|) on basis functions into which the large component of the one- 
fermion 4-spinor is expanded. Hence, the model spaces for the large and the small 
components are generated in terms of this transformation. The advantage of this 
form of the kinetic-balance condition is that the large-component and small- 
component model spaces remain normalized, ft is also possible to transform the 
Dirac Hamiltonian and then form identical model spaces for the large and small 
components. The transformed Dirac Hamiltonian is the so-called modified Dirac 
Hamiltonian [46] and is the basis of orbital-based exact-decoupling methods [;8|. 


The kinetic-balance condition in Eq. ()5]) also ensures the correct non-relativistic 
(NR) limit for c —» oo. The Rayleigh quotient of Eq. (]4]) yields in the limit 
c —> oo the non-relativistic Schrodinger energy: 


E] 


NR 


= lim 

c—yoo 


(V ? ( r> ) I hp - (3mc 2 | j}{r)) 
(il>(r) | i/>(r)) 


p 

2m 


+ R 


l l 


(13) 


where |Z) denotes the (scalar part of the) large component of the spinor after 
taking the limit. 
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2.2 Many-Fermion Dirac Hamiltonian 


The relativistic first-quantized many-fermion Hamiltonian (with positive-energy 
projection not explicitly shown for the sake of brevity) reads, 

N 

= Y, h ^) + w{N) ( 14 ) 

i =1 

with 


h D (i) = 1 4 (1) <g) • • • ® U(i - 1) <g) h D <g) 1 4 (i + 1) <g) • • • ® l 4 (iV) , (15) 

where /i£> is the one-fermion Dirac Hamiltonian of Eq. © taken for fermion i 
and VF (jV ' ) describes the interaction of all pairs of the N fermions. The wave 
function for N non-interacting fermions, i.e., ) = 0, can be constructed as 

the direct product of one-fermion 4-spinors r ip i ( y r i ), 


y{r) = VhOi) ® ® ^i(ri) <g>... <8> ^ N {r N ) , (16) 


which can be antisymmetrized to fulfill the Pauli principle. Now, r = (rq,..., rjy) T 
collects all N one-fermion coordinates. In the case of two fermions, we have the 
direct product of two basis states 


VhOi) <g> OV 2 ) 






VlVO^M 

VlVOVlVO 

VIVO 


^(r 2 ) 


0iV00!V2) 

0i 2 ( r i) 


VlVa) 


VlVO^fo) 

V’f 1 (n) 

<g) 

V’fVa) 


VlVOVlVa) 

VfVO. 


_0!V 2 )_ 


vfVO^fo) 


VfVO^fM 


(17) 


The superscripts T and V indicate large and small 2-spinors, respectively, as 
before. The number attached to these letters indicates the element of a 2-spinor. 
For instance, the elements of the large-component 2-spinor are denoted as 


WO 


VIVO 

vivo 


(18) 


A basis-set expansion of an 77-fermion wave function may be constructed to be 
consistent with the model space 


n {N) _ 0 . . . 0 0 ... 0 n s - s , 


(19) 
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( 20 ) 


where each 'ZZ Al ' Ajv is constructed from the one-fermion model spaces, 

n Xl - XN = n Xi ® ...®'h X n , 


with Ai,..., Ajv G {l, s}. The highlighted spinor components in Eq. (1T7|) are those 
contained within the model space 'H 11 . We recognize that the wave function in Eq. 
(fT6l) and the model space in Eq. (fT9]i are not compatible since it is not possible 
to partition Eq. (ITU]) in terms of the one-fermion model spaces. However, we can 
reorder the spinor elements of the wave function as 

P T 'f’{r) = i/hOi) El... B B B i/) N (r N ) (21) 


where IE is the Tracy-Singh product and P is a permutation matrix (see appendix 
IA.1I for further details). Then, our two-spinor example reads 


V>i(n) B^ 2 (r 2 ) = 


il>[(ri) ® il> l 2 (r 2 ) 

0 ^ a 2 {r 2 ) 
V’iCn) ® 1> l 2 (r 2 ) 
*l> a i(ri) ® ^ s 2 (r 2 ) 


^ 2 (r*i)4 1 (r* 2 ) 

Vf(n)V4V) 

*l>[ 1 (r 1 )i/j?{r 2 ) 

(ri)^| 2 (r 2 ) 


( 22 ) 


The spinor components highlighted in Eq. (12 2 p are those contained within the 
T-L 11 model space as in Eq. m We see that the wave function in Eq. Q 22]) can be 
partitioned such that the individual components are part of the different model 
spaces in Eq. (fT9]h 

[P r T f (r)j Ai '" Aa = <g>... ® ^(ri) ® ® ^(r N ) , (23) 

where Xi,..., X N e {Z, s} as in Eq. (1201) and antisymmetrization will be required. 
The Hamiltonian is transformed accordingly (cf. Eq. (IPPj) in the appendix) 

N N 

H [ ^ s = P T H ( ^ ) P = J2 P T h D {i)P + P T W {N) P = Y h DTS (i ) + P r W {N) P 

i=l i =1 

(24) 

with 

hoTsix) = hoijXjP = 14 ( 1 ) E3 ■ • ■ E l 4 (i — 1) E hp E 1 4 (i -1- 1) E • • ■ E l 4 (iV). 

(25) 
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The potential-energy operator W l - N ' 1 will be invariant under this transformation 
if only the instantaneous Coulomb interaction is considered as it is a diagonal 
matrix with identical entries. The situation is more complicated when magnetic 
interactions are taken into account. An 7V-fermion wave function for 1/2-fermions 
can then be partitioned in terms of the model space into 2 N components each of 
dimension 2 N , 





" \ L 

■0 ' 

T'(r) = 


= 

|Ai.. 

• Ajv) 


V s - s {r) 


\s . 

. .s) _ 


Note that a related reordering of the Hamiltonian similar to Eq. (I24p is key for the 
quaternion formulation of four-component self-consistent field algorithms j.4B|- 


3 Exact Two-Particle Kinetic-Balance Condi¬ 
tion 


In this section, we derive the kinetic-balance condition for explicitly correlated 
basis functions for a system of two fermions. According to Eq. (HUP the model 
space takes the form 

n {2) = n tl © u ls ® n al © n ss , ( 27 ) 

where the four subspaces are formed from the single-fermion model spaces 1-L 1 
and 'H s : 


u 11 = n l ® n l , 

(28) 

u ls = u l ® 

(29) 

n sl = h s ®u\ 

(30) 

u ss = u s ®u s . 

(31) 


Each model space in Eqs. (j28]) - (j3T]) is assigned to one of four components in 
the 16-component wave function. The structure of the Dirac Hamiltonian has to 
respect the structure of the Tracy-Singh product (see Eq. (197|) in the appendix) 
to match the partitioning of the wave function according to Eq. (HMlh We then 
obtain the following block structure for the two-fermion Hamiltonian defined in 






Eq. m- 



c{<r 2 ' P 2 ) 

V + W — 2mi 2 c 2 l4 
(32) 


c(crf } ■ Pl ) 


where we introduced the four-dimensional unit matrix I 4 to highlight the di¬ 
mension and V = [V(t*i) + E(r 2 )]l _4 to yield a four-dimensional respresentation 
of the external potential-energy operator. Moreover, we assume that V and 
also the four-dimensional fermion-fermion interaction operator W are diagonal, 
which does not hold if magnetic and retardation effects are considered for the 
interaction of the two fermions (hence, we apply the compact notation 'W 1 for a 
4x4 matrix operator describing the Coulomb interaction of two fermions only). 
If this assumption is not made, rather complicated expressions will emerge for 
a magnetically balanced, explicitly correlated basis. In particular, the zero en¬ 
tries in Eq. (152|) that represent the cases with a large and small component in 
the bracket per fermion would carry the magnetic fermion-fermion interaction 
(as expressed, for instance, in the Gaunt or Breit operators). As we will later 
make an assumption that all potential energy contributions are small compared 
to the rest energies of the fermions, we aim at a kinetic balance condition free 
of any reference to a potential energy operator in analogy to the orbital-based 
two-fermion case. 

Note that we have also introduced an energy shift of the whole spectrum in Eq. 
(f32ll by — mi 2 c 2 with m i2 = mi + m 2 . Moreover, we absorbed the direct products 
into ct { - 2> as 



(er x <S> 1 2 , cr y ® 1 2 , cr z ® 1 2 ) T , 


(33) 


and 



(I 2 ® cr x , 1 2 <g) cr y , 1 2 <g> ct 2 ) t . 


(34) 


The idea of kinetic balance is to relate the small-component one-fermion model 
spaces to their large-component one-fermion model spaces in the eigenvalue prob¬ 
lem 



( 35 ) 
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This leads to a system of four equations, analogously to Eqs. (0 and (0, 

0 = (V + W - E1 4 )V 11 + c(of } ■ p 2 )^ ls + c(of } ■ pjT' 8 ', (36) 

0 = c(of } ■ p 2 )V u + (V + W - EU - 2m 2 c 2 l 4 )^ 8 + c(<r[ 2) • Pl )^ 88 , (37) 

0 = c(<r| 2) • Pl )¥“ + (V + TE - £1 4 - 2mic 2 l 4 )^ 8 ' + c(<t), 2) • p 2 )^ 88 , (38) 

0 = c(.tS 2) ■ p 4 )^ ls + c(crf } ■ p 2 )V sl + (V + W — El 4 - 2m 12 c 2 U)^ ss , (39) 

where we have suppressed the coordinate dependence of the 4-spinors and will 
continue to do so where convenient. We eliminate one of these four equations 
because we search for a relation between the four four-dimensional components of 
the wave function which we can then apply as a constraint on explicitly correlated 
basis functions. As in the case of a single fermion, we eliminate the energy E 
from the equations by approximating 

[2m,c 2 + E\ 1 4 — V — W 2m*c 2 l 4 (40) 

where rrii G (mi,m 2 ,mi + m 2 }. Similarly to the one-fermion case, Eq. (fTUj) . we 
assume that this approximation remains valid and a variationally stable many- 
particle basis set can be derived. 

We eliminate the first equation, Eq. fl36]h from the system of equations since it is 
the only equation where 2mjC 2 does not occur so that Eq. (140]1 cannot be applied. 
After applying Eq. (140]) to Eqs. (1371) (T39|l . we find the following relations among 
the four components of the two-fermion wave function: 

0 w c((t ( 2 2) ■ p 2 )^ u — 2?7l 2 c 2 ^ s + C (o-f ) -p 1 )^ 88 , (41) 

0 » c(of } ■ Pl )V u — 2mic 2 ^ + c(4 2) -p 2 )^ 8S , (42) 

0 « c(er^ 2) ■ Pl )\I> Zs + c(<t 2 2) • p 2 )\I> 8 * — 2mi 2 c 2 \E f88 . (43) 

The matrix form of this under-determined system of linear equations can be 

interpreted as the augmented form of a linear system with a unique solution: 



( CT 2 2) ' Pi) 

—2m 2 cl 4 

o 4 

-(of 1 - Pi) 

{1} 

A = 

( ct i 2) ' Pi) 

o 4 

—2?n 1 cl 4 

“(°2 2) ‘Pi) 

{2} 


°4 

( ct i 2) ’Pi) 

(o- 2 2) • p 2 ) 

2mi 2 cl 4 

{3} 


The augmented form of linear systems and row reduction are explained in some¬ 
what more detail in appendix IA.21 The number in curly brackets on the right- 
hand side counts every row. ft will be used to express the manipulations in the 
row reduction below. 
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There is no row-reduced echelon form for the augmented form in Eq. (1441) . The 
lack of a multiplicative inverse of the differential operator prohibits setting the 
leading element of each row of the row-reduced echelon form to 1 (see Eq. (11031) 
in the appendix) and therefore to relate A> u (ri,r 2 ) 1 dd s (ri,r 2 ), and '& sl (ri,r 2 ) 
to T ss (r 1 , r 2 ). However, we are able to find a similar form with pairwise relations 
between d' ss (r 1 , r 2 ) and the other three components. These individual steps are 
to be taken in order to obtain this modified row-reduced echelon form: 


1. Insert ((t (2) ■ ^{l} - (cr(, 2) ■ p 2 ){2) into {2}: 


to 

—2 m2cl 4 

o 4 

-(°i ] -Pi) 

0 4 

— 2m 2 c(<Ti ■ p x ) 2mic(cr 2 ■ p 2 ) 

[pi - pi] 14 

0 4 

(^1 -Pi) 

(.*?’ ■ P 2 ) 

2mi 2 cl 4 


{ 1 } 

{ 2 } 

{3} 


2. Insert {2} + 2 ttt,2c{3 } into {3}: 


(aqd ■ p 2 ) —2m 2 cl± 0 4 

1-1 s 

_1 

0 4 —2m 2 c(cr^ ■ p x ) 2 mic(cr 2 2) • p 2 ) 

[pi - pi] 14 

0 4 0 4 2mi 2 c(<r.f ) ■ p 2 ) 

pi~ pi + 4m 2 mi2C 2 l 4 


{ 1 } 

{ 2 } 

{3} 


3. Insert-— {2} H--{3} into {2}: 

m 2 m 2 


(crjjd ■ p 2 ) — 2m 2 cl 4 0 4 

-Pi) 

0 4 2mi 2 c(crf } • p x ) 0 4 

[pi - Pi + 4mimi2C 2 ] 1 4 

0 4 0 4 2m 12 c(cr < £' > • p 2 ) 

[pi — pi + Am 2 mi 2 c 2 } 1 4 


4. Insert (cr i \ 2> ■ p 1 )mi 2 {l} + m 2 {2} into {1}: 

(2) . „ V~(2) 


m 12 ((T ( 2 ■p 2 ){v{ -Pi) 

0 4 

0 4 


0 4 


2m 12 c(cr l 

0 4 


( 2 ) 


0 4 

• Pi) 0 4 

2m 12 c(cr 2 ) • p 2 ) 


-m 1 p\ 


\pi 


{i} 

{ 2 } 

{3} 


■ m 2 p 2 + 4mim 2 mi2C 2 ] 1 4 
- p\ + 4mimi 2 c 2 ] 1 4 


{ 1 } 

{ 2 } 

{3} 


We arrive at a set of simple pairwise relations between T ,ss (r 1 , r 2 ) and the other 
three components 

-(erf 21 ' Pi)(^ 2 2) ‘ f> 2 ) m i 2 'S' 11 = (mipl + m 2 p 2 2 - Amim 2 m 12 c 2 ) d> ss , (45) 

-2c(crp ) ■ p 1 )?n 12 T ,Zs = (pi - pf - 4m 1 mi 2 c 2 ) T' ss , (46) 

-2c(ct2 2) • p 2 )mi 2 '$! sl = (pi - pi - Am 2 m 12 c 2 ) d> ss . (47) 
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Forming the least common multiple from the operators on the left-hand sides of 
the equations, we can introduce a four-dimensional spinor @(ri,r 2 ) related to 
the ^ ,ss (ri,r 2 ) component, 

^ ss (r 1} r 2 ) = —2cmi 2 (cr ( L 2) • p x )(c t^ ] ■ p 2 )0(n, r 2 ) , (48) 


insert it into Eqs. (H5]) - (jTT)) and eliminate identical terms on both sides. Instead 
of relating the upper component to the lower component, we relate all four four- 
dimensional components of the 16-spinor, 


W)' 


>>' 

V ls (r) 


\ls) 

V sl (r) 


\sl) 

V ss \r) 


ss) 


to a common spinor ©(rq, r 2 ): 


(49) 


| ll) = (mipl + m 2 Pi ~ 4mim 2 mi 2 c 2 ) l 4 ©(rq, r 2 ) = U^ ) @(r 1 , r 2 ), (50) 

\ls) = (T \ C P2 (P 2 - Pi - 4mimi 2 c 2 ) 0(ri, r 2 ) = t/ ; (2) ©(r,, r 2 ), (51) 

/ (2) . \ 

|sZ) = — 1 (p 2 - pi - 4m 2 m 12 c 2 ) ©(rq, r 2 ) = Z7 s ( , 2) 0(n, r 2 ), (52) 

| ss ) = -m 12 (crS 2) ■ £> 1 )(<t), 2) ■ p 2 )0(rq, r 2 ) = r 2 ). (53) 

Here, we have introduced the short-hand notation u[f\ U^\ U^\ and Usf for 
the transformation to kinetically balanced components in analogy to the one- 
fermion case in Eq. (fT2|) . In a subsequent section, we refer to the i-th term in 
the prefactor of such expressions as d[ N ' > with N— 2 for the two-fermion case; e.g., 
d ^' 1 for | si) is then — cr[ 2) •p 1 /2c x 4m 2 mi 2 c 2 . The physical role of ©(rq, r 2 ) will 
become clear when we study the non-relativistic limit (see below). We emphasize 
that ©(t*i, r* 2 ) is in general an explicitly correlated geminal rather than a simple 
orbital product. 


Because of the derivation in Eq. (l48]h T fSS (ri, r 2 ) is uniquely defined by @(r*i, r 2 ) 
up to a constant, i.e., the constant of integration. For square-integrable functions, 
this constant is zero. Hence, cancellation of differential operators is not a problem 
and all components are uniquely determined by ©(r 1 ; r 2 ). 

Finally, we consider fermion exchange symmetry (Pauli principle) for the two 
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identical fermions leading to the relations [38] 

^ H (fT,** 2 ) = -^ H (r 2 ,ri), (54) 

'f' ls ( y r 1 ,r 2 ) = ~'& sl {r 2l r 1 ), (55) 

^ M (fT,r 2 ) = -^ ss (r 2 ,ri), (56) 

which have to be fulfilled in addition to the relations in Eqs. fl50|) - (l53j) . ©( 7 * 1 , r 2 ) 
is antisymmetrized before the components are constructed according to Eqs. 
§UD-§3D because the operators ( cr ^ ■ p-J and (cr^ ■ p 2 ) do not commute with 
the permutation operator which exchanges fermions 1 and 2. 


4 The Non-Relativistic Limit 


The one-fermion kinetic-balance condition yields the correct non-relativistic limit 
for c —y 00 . This is a key requirement ensuring variational stability. We therefore 
require any kinetic-balance condition for more than one fermion to yield the 
correct non-relativistic limit. 


Finding the non-relativistic limit for the one-fermion case is fairly trivial. For 
the two-fermion kinetic-balance condition, this is somewhat more involved. In 
order to find the correct limit, we rely on de l’Hospital’s rule for limits, 


Km hd = lim hW, 
x^ry g(x) x->y g'(x) 


(57) 


where f(x) and g'(x) are the derivatives of f(x) and g(x) with respect to x, 
whereas y is the limiting value of x. 


The non-relativistic limit of the two-fermion total energy for a wave function 
kinetically balanced according to Eqs. (j50j) - (l53D , can be taken as a limiting case 
of the Rayleigh quotient 


£ nr = lim 


(V\H 


( 2 ) 

DTS 


l*> 


1 ^) 


(58) 
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For the one-electron part in (\I/| H DTS |\I>) we have 
2 

<^| h DTs(ri ) |^} = (ll\ c{ct ( 2 2) ■ p 2 ) | Is) + (ll\ c(cr[ 2) • p x ) | si) 

1=1 

+ (ls\ c(cr ^ • p 2 ) KO — (^ s l 2m 2 c 2 |Zs) 

+ (Zs| c(<t ( l 2) ■ p x ) |ss) + (sl\ c(cr[ 2) ■ p x ) | ll) 

— {sl\ 2miC 2 |s/) + (sl\ c(cr^ ■ p 2 ) |ss) 

+ (ss| c(cr^ 21 • p x ) | Is) + (ss| c(cr ^ • p 2 ) |sZ) 

— (ss| 2mi 2 c 2 |ss) + (thl V ® 1 4 |\I/) , (59) 

where we have not resolved the potential-energy expectation value for conve¬ 
nience. It must now be noted that 

(ls\ c(cr 2 2 ^ • p 2 ) | ll) — (ls\ 2m 2 c 2 |/s) + (ls\ c(cr[ 2) ■ p 2 ) |ss) = 0, (60) 


which can be shown by exploiting Eqs. (115]) and (1461) to replace \ll) and | Is) by 
expressions for |ss). Analogously, we can exploit Eqs. (TTol) (1171) to show 

(sl\ c(cr[ 2 ' ) • p x ) | ll) — (s/| 2m!C 2 | si) + (sl\ c(cr 2 2 ^ ■ p 2 ) |ss) = 0 , (61) 

(ss| c(cr^ . p 1 ) | Is) + (ss| c(<t 2 2 ' ) ■ p 2 ) |sZ) — (ss| 2mi 2 c 2 |ss) = 0. (62) 

Hence, we find for the full Hamiltonian with interacting fermions 

<*l H'Sts I*) = <H| c(of > ' P 2 ) M + W <:(*? ■ Pi) M) + <*| (V + W) ® 1 4 |*). 

(63) 


We now apply de FHospital’s rule to Eq. (158|) by taking the fourth-order derivative 
with respect to c of both the numerator and the denominator: 


E nr = lim 

C—> OO 


dd_ 

dc A 


c(cr 


S 2) ■ Pi) | Is) + (U\ c(«r< 2 > ■ p 2 ) I si) + <*| (V + W) ® 1 4 |*>} 


I#) 


= lim 


(@| 192m 2 2 m 1 m|p 2 l 4 + 192m 2 2 m 2 m 2 p 2 l 4 + 384m^ 2 m^m 2 (y + W) + 0(c 2 ) |@) 


384mf 2 mfm| (© |@) 


(64) 


The potential energy term, V + W, may also contain contributions depending on 
c, but these contributions are of second or higher order in c^ 1 . When taking the 
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limit, they are all zero and we find the limit to be a simplified Rayleigh quotient 
depending on @(ri, r 2 ) 


-^nr — 


/© 

f +V + W 

©\ 

\ 

2mi 2 m 2 

/ 


(0i©> 


(65) 


where V and W are the limiting values with c —> oo for V and W, respectively. 
In Eq. fl65[) . we obtain the Schrodinger energy and therefore the correct non- 
relativistic limit. The limit also identifies the four-dimensional spinor ©(7*1, 7*2) 
as the non-relativistic two-fermion Schrodinger wave function (note that this 
function still features a four-dimensional spinor structure as it accounts for the 
spin of two electrons). 


It is interesting to note that the value of the non-relativistic, c —>■ 00 , limit is 
determined by the leading terms in c of the three components \ll ), | Is), and | si) 
in Eqs. (I50ll - fl52|l define the non-relativistic limit when we apply de l’Hospital’s 
rule. These leading terms are 


\ll) (c 2 ) : 

-4m 1 m 2 m 12 c 2 ©, 

(66) 

1 Is) (c) : 

-2mim 12 c((T {2) ■ p 2 )@, 

(67) 

and 



1 si) (c) : 

-2m 2 m 12 c{crf ) ■ p x ) ©. 

(68) 

We also note that Eqs. (I66l)-fl68l) 

are related to Eq. ((SJ). If we 

apply Eq. dHJ for 


particles 1 and 2 subsequently to ©(r) and then multiply by 4mim 2 mi 2 c 2 , 


I H) (c 2 ) 

I Is) (c) 

I si) (c) 
\ss) (1) 


—> 4mim 2 mi 2 c 2 


© 

(°2 2) • P2) Q 
2 ?n 2 c 

(°~i 2) -Pi) e 
2 m ±c 

(^l 2 ’ -Pl){<*2 2) -P2) q 

4'm!m 2 c 2 


(69) 


we obtain the expressions of Eqs. f!66p — (I68p . Hence, we have shown that the 
one-fermion kinetic-balance condition in Eq. (E]) is sufficient for obtaining the 
correct non-relativistic limit for a two-fermion system. At first sight, this seems 
reassuring as obtaining the correct non-relativistic limit has been connected to 
variational stability for orbital-based theories (see, e.g., Ref. M). However, the 
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one-fermion kinetic-balance condition may not be sufficient to ensure variational 
stability in the general case considered here dHHSHE]. Accordingly, the non- 
relativistic limit will then not be a sufficient, albeit a necessary condition for 
variational stability. 


5 Kinetic-Balance Condition for More Than Two 
Fermions 


The derivation presented in section [3] can also be applied to systems of more 
than two fermions, and thus establishes in its full form an exact kinetic-balance 
condition for general (non-separable) iV-particle basis functions. How such a 
generalization could be achieved for the approach of Pestka and co-workers [1311361 
07] is not obvious and was not discussed in their papers. In our ansatz, we obtain 
rather lengthy expressions for three fermions, which we refrain from presenting 
explicitly for the sake of brevity. The resulting expressions can, however, be 

/o\ 

expanded into a polynomial with respect to c. The individual terms d\ ' (c) of the 
prefactor of the 3-fermion 8 -spinor ©(rq, r 2 , 7 * 3 ) feature the important property 

df\c) = f(mi,ffl 2 ,m 3 ) x c (#-«-"-«0 (cr^ 3) ■ Pi) u (ctP ■ P 2 ) v (<rP ' Pa) w (70) 

where we have omitted to indicate that each dP (c) will be different for different 
sectors | III), \ lls), |/ss), and so forth and depend on u,v,w. The positive semi- 
definite exponents u,v,w obey the constraints 0 < (u + v + w) < 7 and we 
have 

= ( cr x ® 1 4 , cr y ® 1 4 , <j z ® 1 4 ) T , (71) 

crip = (I 2 <E) cr x <g> 1 2 ,1 2 ® cr y ® 1 2 ,12 ® cr z ® 1 2 ) T , (72) 

crip = (1 4 ® cr x , 1 4 < 8 ) cr y , 1 4 ® cr z ) T . (73) 

The multiplicative prefactors k\ (mi,m 2 ,m 3 ) depend on the masses of the in¬ 
dividual fermions and the kinetic-balance conditions simplify significantly if all 
three fermions have equal masses. 

Eq. (1701) shows that the explicitly correlated kinetic-balance condition for three 
particles contains the momentum operator to the power of seven, which is unfa¬ 
vorable from a computational point of view. However, we can observe that the 
power of the momentum operators decreases with increasing orders of c. The 
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leading terms with respect to c are the one-fermion kinetic-balance terms and 
ensure the non-relativistic limit. 

For the assessment of the general properties of an iV-fermion kinetic-balance 
condition, let us first re-write the two-fermion kinetic balance condition, Eqs. 
(I50j) -(j53 ]) . in a general form similar to Eq. ([TO]) : 

d?\c) = x c< 4 -’-’V 1 2 > • Pi r(a?> ■ p 2 )« , (74) 

where the multiplicative prefactors m 2 ) depend on the masses of the two 

fermions and the positive semi-definite exponents, u and v, obey the constraints 
0 < {u + v) < 3. 

By comparing the results for two- and three-fermion systems, Eqs. (1701) and (1731) . 
we obtain for the IV-fermion case: 

N 

df>(c) = ..., m N ) x c {2N ~ u) JJ(^- Ar) • Pj) Uj (75) 

3 = 1 

where we skipped the explicit derivation. The power of c, 2 N — u, and the power 
of the c t ^ ■ Pj operator, Uj, are determined in the exact kinetic-balance solution 
by 

N 

o < u = < 2N + 1- (76) 

3 = 1 

High powers of the momentum operator is unfortunate from a computational 
point of view, but with the complete set of kinetic-balance conditions at hand for 
any set of non-separable iV-particle basis functions, Eqs. d75|) and f!76]) . one may 
introduce a hierarchy of approximate kinetic-balance conditions and investigate 
their properties systematically. 

As an example, we present the approximate kinetic-balance condition for a three- 
electron system (in Hartree atomic units and with m e = 1 for the electron mass) 
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where only the leading terms in c are included: 


1 III) 

= (48c 6 l 8 -14((^ 3) - Pl ) 2 + ( 

o'? ' P 2 ) 2 + (o"3 3) • 

p 3 ) 2 ) c 

■ 4 J © 

77 

$ 

s? ~ 

77 

III 

77 











(77) 




| lls) 

= ( ct 3 3) ’Ps) ((— (^i 35 -Pi? ~ 

■ P 2 ) 2 - 7(of } 

• P 3 ) 2 : 

)c 3 + 

24c 5 V 

) ®M = 

rj(3) 

U Us 

©( 









(78) 




| Isl) 

= ( ct 2 3) -P 2 ) ((-(^i 3) 'Pi) 2 ~ ‘ 

7(o-2 3) -P 2 ) 2 - ( ct 3 3) 

• P 3 ) 2 : 

)c 3 + 

24c 5 1 8 N 

) ®M = 

f/( 3 ) 
U lsl 

©( 

>)» 








(79) 




/ss) 

= 12 (c4 3) •p 2 ){tr■ 

P:s)c 4 ®(r 

) = l/<2e(r), 




(80) 




| sll) 

= ( ct i 3) ' Pi) ((— 7(cr 

! 3) -pi) 2 - 

/_(3) _ \2 / (3) 

( cr 2 ‘ P 2 ) (.°3 

• P 3 ) 2 : 

)c 3 + 

24c 5 

) ®M = 

f/( 3 ) 
WZZ 

©( 

>)» 








(81) 




s/s) 

= 12(<tS 3) •p 1 )(<r^ 3) • 

p 3 )c 4 G(r 

) = u£>e(r), 




(82) 




ss/) 

= 12(crj 3) •p 1 )(o-f ) • 

P 2 )c 4 &(r 

) = U$e(r), 




(83) 




sss) 

= ■ Pl ){(7 { 2 3) 

■P 2 )(°'3 3) ' 

Pi )c 3 0(r) = USl&(r), 



(84) 




with the cr- 3 ' 1 ■ p t operators dehned in Eqs. (1711)-(1721). 

0(r) 

with 

r = (r 1 ,r 2 ,r 3 ) T 





is the non-relativistic limit of ^(r). We see that the lowest order of c to consider 
is 3 due to the |sss) component. Eqs. (177 j) - ([M|) can be considered as a minimal 
explicitly correlated kinetic-balance condition for a three-electron system. 


6 Basis-Set Expansion and Numerical Results 


In practice, a many-particle wave function can be expanded into a basis set 

*w=EE c ^V) (85) 

i AG A. 

where Cf are the expansion coefficients and 3?^(r) are the basis functions. A 
is the set of all component-index strings consisting of Vs and s’s according to 
Eq. (12611 . i.e., it is the set of 2 N strings of such indices of length N for an N- 
fermion basis function. For the sake of clarity, we explicitly provide the basis 
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functions for the two-fermion case, 


rsjn 


" 0 ' 


' 0 ' 


’ 0 ' 

0 

, = 

z ls 

, &? = 

0 

~ sl 

, and = 

0 

0 


0 




0 

_ 0 _ 


0 


_ 0 _ 


L^: J 


•— [( ~ [g ~ gl ~ gg 

where the four-dimensional basis functions , <&■ , <&• , and are promoted 
to 16-dimensional functions for a compact notation of the expansion in Eq. (l85lh 
note that we write ’O’ in Eq. (l86|i to indicate four-dimensional null vectors for 
the sake of brevity. Eventually, these four-dimensional basis functions are to 
be expressed in terms of basis functions ©j that represent the common non- 
relativistic limit © according to the analysis presented above. 


A transformation, similar to that in Eq. (TT2T) for the one-fermion case, can be for¬ 
mulated for the explicitly correlated kinetic-balance condition in the two-fermion 
case, 


[/( 2 ) — 
u KB — 


u, 


( 2 ) 


id, 21 1 


0 4 


0 4 


O 4 


0 d 


U, 


( 2 ) 


Is 


id, 2 ’1 


Oj 


04 


04 


0 4 


U. 


( 2 ) 


sl 


id, 2) i 


0, 


0 4 


0 4 


0 4 

TjW 

u ss 


( 2 ), 

ss 


|K 


(87) 


in the notation introduced in Eqs. (I50j) - (l53]) and with a normalization introduced 
for each basis-function component according to 


K« 2) | ^ ■ K^IQi) 




4 2M 


7 ( 2 ) 1 


( 88 ) 


and so forth for the other A; note that we dropped the basis-function index on the 
left-hand side for the sake of brevity. Essentially, we normalize each component 
of each basis function individually to ensure numerical stability when solving 
the eigenvalue problem. This procedure can be understood as the relativistic 
counterpart of the quasi-normalization in pre-Born-Oppenheimer theory [5D] . 
Hence, explicit normalization of a trial wave function has to be taken into account 
when the energy is calculated. 
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/ O \ 

In full analogy to the two-fermion case, we construct U y K 3 from Eqs. (177l) - fl84l) . 
In general, the iV-fermion trial wave function is expressed in terms of the trans¬ 
formation as 


V(r) = 


E 


ttW 

r<ll...l u iu.,.1 („'i 

\ u u...i\ 


c? 


u 


(N) 


K 


-0.fr) 


fiss...s u ss...s f r ) 

* leti ,( ) 


E C ‘- Pkb- Ml , (89) 


where the /lU^ 


are the entries of the diagonal matrix U 


(TV) 

KB 


normalized by 


|t/f>1 - (90) 

(with the index i dropped for the sake of brevity as before). The vector con¬ 
tains the non-relativistic limit, ©., 2 N times as entry, i.e., 0 ,- A ^ = (©,, ©.,..., ©,.). 


6.1 Numerical Results 


As an example, we present numerical results for a standard two-electron system: 
two electrons moving in the central potential of a helium nucleus within the 
Born-Oppenheimer approximation. Our starting point is a non-relativistic basis 
set, which corresponds to L = 0 total spatial angular momentum, p — +1 parity, 
and S = 0 total electron spin quantum numbers, and which is antisymmetrizd 
according to the Pauli principle: 


©*(*•) 


©i(n,r 2 ) 


1 

7! 


i 

0 


0 

1 



where we inroduced the explicit form of the spin functions 



a = 



and (3 


0 

1 


(91) 


(92) 


In this notation, the four-dimensional structure of the non-relativistic limit is 
highlighted in agreement with Eq. (fSSD - In Eq. (1911) . the spatial part can be 


20 











any non-separable two-particle function and in our calculations it is an explicitly 
correlated Gaussian function with L = 0 and p = +1, 

®/0T,r* 2 ) =exp ^-^r T (A ?; ® 1 3 ) , (93) 

where r = (rq,r 2 ) T and the elements of the symmetric, positive definite matrix, 
A,- l G M 2x2 , are parametrized by 

{Ai} kl = S k i exp (a k i,i) + 0.1(4/ - 1) exp (~a k i,i) with k, l G {1, 2}. (94) 

The a k i,i values are optimized stochastically to minimize the relativistic energy. 
Trial values for a k i,i were generated from a normal distribution as in Ref. [.50] (see 
also references therein). The optimized parameter values of Ai are deposited in 
the supplementary information. 


With Eq. (159]) (see also Eqs. (15U]) - (155]) for the two-particle case) we generate a 
kinetically balanced basis set from ©j(r), Eq. (191]) . for the relativistic calculations 
and minimize the Rayleigh quotient, Eq. (J3J), 




-ef(r) 

ix( 2 ) 

DTS 

GAefV)) 

Ey c-c s (t/£> ■ ©f> 

(r) 



(95) 


with respect to the expansion coefficients Cf by solving the generalized eigenvalue 
problem 


HC = ESC . (96) 

In Eq. (199]) . the Hamiltonian matrix, H, has a block structure with H,j = 

(■ i,j = 1, 2 ,..., n) and simi- 


v IS ■ ©fV) 


H 


( 2 ) 

DTS 


u\& ■ ef (r) 


X2 J 


larly the overlap matrix, S, contains S t] = (U^ ■ ©) '(r) 

M 2 ‘ v x2 jv j _ 1 ^ 2 , ... ,n) for n basis functions and for two electrons, N= 2 . Ac¬ 
cordingly, C G ]R r!2iVxn2jv is a matrix containing the expansions coefficients C( A 
and E is an ( n2 N ) - dimensional diagonal matrix with the energies on its diagonal. 


v(2)/ 


r/( 2 ) . © 

u KB 


( 2 ), 


r) 


The ground-state energy eigenvalue of the helium atom is obtained from Eq. (199]) 
by direct solution of the generalized eigenvalue problem in the stochastically 
optimized basis set (see Table [I]). The non-relativistic energies, also given in 
Table |H were obtained from the generalized eigenvalue problem solved for the 
Schrodinger Hamiltonian in the basis of the non-relativistic basis functions of 
Eq. (191]) . containing the parameters obtained in the relativistic calculations (see 
the supporting information for details). As it can be seen from the data in Table 
HI both the relativistic and the non-relativistic energies converge with increasing 
basis-set size towards the reference data in a variationally stable fashion. 
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Table I: Ground-state energy of the two-electron helium atom with fixed nucleus obtained for 
increasing basis-set sizes. A E R and AEnr are the differences between the calculated relativistic 
and non-relativistic energies, respectively, with respect to the reference values, n is the number 
of basis functions, 0*, defined in Eqs. m and (1931) . The parameters of the basis functions are 
deposited in the Supplementary Material and the value for the speed of light was set to 137.0359895 
atomic units. 


n 

E r [E h ] 

A£ r [E h ] 

-Enr [Eh] 

A-Enr [Eh] 

10 

-2.89757665 

0.00628019 

-2.89744422 

0.00628016 

20 

-2.90288205 

0.00097479 

-2.90275061 

0.00097377 

50 

-2.90382266 

0.00003418 

-2.90369103 

0.00003335 

100 

-2.90384822 

0.00000862 

-2.90372140 

0.00000298 

200 

-2.90385566 

0.00000118 

-2.90372429 

0.00000009 

300 

-2.90385674 

0.00000010 

-2.90372430 

0.00000008 

Ref. [36] 

-2.90385684 

Ref. [oT 

-2.90372438 



7 Conclusions 


The kinetic-balance condition for the one-fermion case ensures variational sta¬ 
bility in orbital-based approaches to first-quantized relativistic many-fermion 
theory. In the present work, we derived a kinetic-balance condition for gen¬ 
eral, non-separablc IV-particle basis functions. Similarly to the derivation of a 
one-particle kinetic-balance condition, we set out from the assumption that the 
potential energy contributions are small compared to the rest energies of the 
fermions. We arrived at an IV-particle kinetic balance condition by combining 
the well-known multiplication properties of the Pauli matrices with the row- 
elimination approach of solving linear systems of equations. In agreement with 
the one-fermion case, the IV-particle kinetic-balance condition also ensures that 
the correct non-relativistic limit is obtained for an infinite speed of light. It had 
been anticipated, however, that the IV-particle kinetic-balance condition provides 
better stability when solving the first-quantized Dirac Hamiltonian variationally 
with an explicitly correlated basis set, and hence suggested that the requirement 
of matching the non-relativistic limit is a necessary but not a sufficient condition. 

We demonstrated that the variational solution of the Dirac equation is stable for 
the ground state of the two-fermion helium atom when a relativistic basis set 
is generated from explictily correlated Gaussian functions using the IV-particle 
kinetic balance condition for N= 2. 
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Concerning the general applicability of our results, the theoretical expressions 
and our preliminary investigations show that the direct use of the full TV-particle 
kinetic-balance condition becomes tedious and computationally expensive for 
more than two fermions. However, it might be possible to reduce the com¬ 
putational cost by systemtically eliminating terms of high order in momentum 
operators from the exact expressions and, at the same time, retain variational 
stability for the solutions. A systematic investigation of the variational stability 
under such approximations is beyond the scope of the present paper and left for 
future work. 
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A Appendix 

A.l Tracy—Singh Product 

The Tracy-Singh product [52j is defined as 


(B n <S> C uv ) 


(B\ n C U v) 


A tsp — BBC — [{Bij <S> C uv )\ 


(97) 


(■®ml ® ^uv) * ’ * (-^mn ® ^uv) 


where B = [B^] and C = [C uv \ are two matrices of dimension (■ mxn ) and ( pxq ), 
respectively. They are partitioned block-wise in terms of the matrices B t j and 
C uv . A tsp is a matrix of dimension (rnp x nq). It is partitioned block-wise with 
the elements being the matrices (B^ <E> C uv ). The Tracy-Singh product may be 
considered a more general form of the Kronecker product 



Akp — B $5 C — [(^pc-uu)] 


(98) 
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where bij and c uv are the matrix elements of B and C, respectively. A^ p is a 
matrix of dimension (rnp x nq). The two matrices A tsp and A^ p are identical in 
the case that B and C are not partitioned (or partitioned into (1 x 1) blocks). 
Generally the two products are related through a permutation of the row and 
column space of either matrix joTH Iq5] 

P 1 \B 1 <g> ... <g> B n )Q = Bl .. MB n (99) 

where P and Q are the permutation matrices for the row and the column space 
and n is the number of matrices involved. For vectors v t , we find the relation 

P T (v i <g)... ® v n ) — Vi IE ... IEI v n . (100) 

The partitioning of the matrices and vectors depends on the permutation matrices 
P and Q. If all matrices are square and symmetrically partitioned, the two 
permutation matrices are identical [53] P = Q and the two products are related 
through a unitary transformation. 


A.2 Row Reduction and Row Reduced Echelon Form 

Systems of linear equations are conveniently solved by first representing them in 
matrix form 


Ax -6 = 0 , ( 101 ) 

where A is a matrix containing the linear factors, a? is a vector and contains 
the values which are to be determined and b is a vector containing the constant 
factors of the linear system. A reliable method of solving such a linear system 
is row reduction, i.e., Gaussian elimination. It involves performing a series of 
operations on the augmented form 

A aug = [A\b] (102) 

until it is in row-reduced echelon form. The row-reduced echelon form is 

A Iie = [l | b'] (103) 

for systems with a unique solution. Possible operations are permutation of two 
rows, multiplication of individual rows with a constant scalar factor and evalu¬ 
ating the difference of two rows. 
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